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1 Introduction 


MARSYAS is a computer-aided control system design package for the simulation and analysis of dynamic 
systems. In the summer of 1991 MARSYAS was updated to allow for the analysis of sampled-data systems 
in terms of frequency response, stability, etc. This update was continued during the summer of 1992 in order 
to extend further MARSYAS commands to the study of sampled-data systems. Further work was done to 
examine the computation of OPEHAT transfer functions, root-locii and w — plane frequency response plots. 


2 Sampled-data systems with feed-forward coefficients 

Consider a sampled data system whose state- space equations are 


x c (t) = A c x c {t) + B cd y d (kT) + B c u(t) (2*1) 

y(t) = C c x c (t) + D cd y d (kT) + D c u(t) (2.2) 

u d {t) = C dc x c (t) + D dd y d {kT) + D dc u{t) (2.3) 

with corresponding discrete time subsystem 

x d (kT + T) = A d x d (kT) + B d u d (kT) (2 A) 

y d (kT ) = C d x d (kT) + D d u d (kT). (2.5) 


(Purely discrete time inputs r d (kT) to the discrete time system may be incorporated into the above equations 
by augmenting the input vector u(f) and the feed-forward matrix D dc .) The feed-forward coefficient D d 
causes the states x c (t) to depend not only on their continuous values through the coefficient matrix A c , but 
also on their sampled values through the coefficient matrix B cd D d C dc * As of MARSYAS version 6.0.5, this 
dependence is dealt with by creating a matrix ZX = B cd D d C dc in the linearized continuous time dynamics. 

Unfortunately, MARSYAS version 6.0.5 does not correctly treat the case where there are feed-forward 
coefficients in both the continuous time and discrete time blocks (i.e. D dd ^ 0 and D d ^ 0). In this case, 
MARSYAS infers an algebraic loop where there is none. For example, consider the sampled-data system 
with continuous time block 

y(t) = y d {kT) u d (i) = u(t)/2 - y d {kT) 

and discrete time block 

yd (kT) = k d U d (kT). 

The code has no (explicit) states, yet because of the direct feed-forward gains (i.e., the U D” matrices) in the 
discrete and continuous time blocks of the linear system, a “hidden state” associated with the A/D converters 
in the system manifests itself at sampling times. The MARSYAS-6. 0.5-generated system equations are 

W [ 1] = ( 2 . 500000E-01 ) * U [ 1] 

while the simulated output values display “spikes” at sampling times and the values between sampling times 
incorrect. There is a specific order in which the simulation updates must occur at each “scheduled” run of the 
discrete time systems: (1) Evaluate the continuous time system states and algebraic outputs. (2) Update the 
discrete time system states and algebraic outputs. (3) the algebraic outputs of the continuous time system to 
reflect the discrete time system changes. Once these three steps are completed, the simulation may proceed 
with the next integration step. Similarly, the analysis of discrete time systems must be modified as follows. 

Lemma 2.6 Let a sampled data system be defined by the equations (2A)-(2.5) with sampling time T given. 
Then 

u d (kT) = C de x c (kT ) + D dc u(kT) + D dd y d (kT - T); 
that is, the discrete-time output y d becomes a system state when D dd ^ 0. 


XXI- 1 




Figure 1: OPEHAT linearized system 


Based on the above operation, the vector jj<j( 0 can be used as a vector state that replaces the “ZX” input 
used in MARSYAS version 6.0.5 as follows. Let x dy {kT) = y d {kT - T ). Then the overall system simplifies 
to 


*«(<) 

Wrf(t) 

x d {kT + T) 
x dy (kT + T ) 

Vd{kT) 


A c x e (t) + B cd y d {kT) + B c u(t) 

C e x c (t ) + D cd y d (kT) + D c u(t) 

C dc x c (t) + D dc u(t) 

A d x d (kT ) + B d u d (kT) + B d D dd x dy {kT) 
C dc x d (kT) + D d D dd x dy (kT) + D d u d (kT ) 
C d x d (kT) -f D d D dd x d y 4- D d u d {kT). 


3 OPEN AT analysis: options and an example 


The OPEHAT command allows MARSYAS users to examine system robustness with respect to an individual 
parameter by breaking a specified signal path in a closed-loop system model and examining the poles/zeros 
of the newly created open loop system. 

The current MARSYAS implementation (and the original MARSYAS implementation) treat the OPEHAT 
command as a special case of the IHOUT command. The manual for the original MARSYAS implementation 
describes an alternate technique for computing the OPEHAT transfer function; this alternate technique forms 
the basis of the root locus calculation, and is thus of interest in both of these calculations. The two techniques 
are described below. A simple numerical example is given that indicates that algebraic loops (i.e., a non-zero 
D matrix in a linearized OPEHAT system) can render the alternate algorithm inaccurate; it is advised that 
the IHOUT approach be used in both OPEHAT and root locus calculations. 


The OPEHAT command can be summarized as follows. Given a closed loop continuous time system 
x = f ( x ’ with algebraic constraint 0 = g(x , u>) , OPENAT k selects a gain value in the system and linearizes 
about an operating point in order to obtain the SISO system shown in Figure 1. The OPEHAT command 
identifies the poles and zeros of the transfer function n(s)/d(s). The IHOUT algorithm computes the poles 
and zeros by computing the matrices (A, B, C, D) that characterized a state-space realization of n(s)/d(s) 


and then computing the (finite) generalized eigenvalues of the matrix pencil 



The alternate closed-loop matrix method is as follows For a fixed linearization, let A(k) be the system 
Jacobian matrix evaluated with K = k\ i.e. A(k) = (A + k/{\ - Dk)BC). Then det(sJ - A(*)) = 
a(fc)(d(s) - kn(s)) where a(fc) -1 is the leading coefficient of (d(s) - *n(s)). Clearly, the poles of the OPEHAT 


system are the eigenvalues of A(0), since det(s7 - A(0)) = a(0)(d(s) - 0 • n(s)) = d(s). 

In order to obtain the zeros of the system, it is required to compute the algebraic gain D and the closed 
loop matrix A( 1) = (A(0) + 1/(1 — D)BC). Let B, C be a column and row vector, respectively, such that 
BC = (D — 1)(.4(1) — .4(0); then the zeros of the OPEHAT transfer function may be obtained by computing 
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INOUT method Closed-loop matrix method 
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Figure 2: Results of tests of openat algorithms 


the finite generalized eigenvalues of the matrix pencil ^ ^ »*]-(;:]) . This is the closed-loop 


-A -B 
C D 

matrix method discussed in the manual for the original MARSYAS implementation. 

The discussion of the closed-loop matrix method in the original MARSYAS manual did not correctly 
treat the case of D ± 0. If D ^ 0, in particular, if D « 1, then the INOUT method is clearly superior to the 
closed-loop matrix method, as shown in the following example. 


Example 3*1 Let the OPENAT linearized system have coefficient matrices A = 
C = [ 2 1 ] , and D — 1 -b e. The exact system zeros can be found to be 


-1 1 
-1 -1 


B 


1 ' 

2 ’ 


A = 



2 

1 + e 



Observe that as e — ► 0, the zeros coincide at A = 3. 

A computed linearization (A, S, C, D) = (A, B } C } D) will typically have roundoff noise in each matrix in 


proportion to the matrix norm; i.e., 


a|| < fi\\A\\ where is a function of machine precision and the 

algorithm condition. The example was selected so that matrix entries were “balanced;” that is, the matrix 
entries were of approximately the same magnitude so that round-off effects on problem condition would be 
reduced. For the purposes of a numerical experiment, small error (« 10” 14 ) was deliberately introduced into 
the problem in order to examine algorithm sensitivity. Since the computed value of Ab is near to the correct 


computed value, this method yields good results. 

The closed-loop matrix method involves computing A(0), A(l), and the open-loop D-matrix. it was not 
assumed that the gain D involved in computing A(F) = A — 1 /(D — 1 )BC is exactly the same as that gained 
from an IMOUT linearization; i.e., INOUT computes (A, £, C, D\) = (A, £, C, £)), while A(l) was computed as 
A(l) = A+ 1/(^2 - 1 )BC. The respective results are shown in Figure 2. It is clear that the INOUT method 
yields superior accuracy; it is important to recognize that the difference in these two methods is a result of 
the difference between 1 /{D\ — 1) and 1/(J?2 — l)j since the estimated values of B and C in the closed loop 
matrix method will be scaled by the ratio (D\ — l)/(D 2 ~ 1) from their “correct” values, it is to be expected 
that the closed-loop matrix method will yield poor results when D « 1. 


□ 


On the basis of the above example and analysis, it is recommended that both OPENAT and root locus 
calculations be performed on the basis of the matrices (A, £,C, D ) obtained from an INOUT linearization of 
the appropriate system. FORTRAN code for the mvzero routine and the associated numerical balancing 
procedure [2] have beed delivered to John Tiller, BOSS, this summer for use in MARSYAS. 
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The calculation of OPEHAT transfer functions and root locus data become much more complicated in a 
sampled-data system. If the OPEHAT operation is performed on a gain k in the discrete time portion of the 
system, then the OPEHAT function may be executed in the normal fashion on the equivalent discrete time 
system. However, if the OPEHAT is performed on a gain k in the continuous time part of the system, then the 
OPEHAT function ceases to have a clear meaning. If an artificial continuous time input it and output y are 
placed around the gain, then there is no rational transfer function from 1 1 to y because of the time-varying 
dependence on the discrete time subsystem. Instead, discrete time input-output pair is inserted about the 
continuous time gain element K. 

This system is not amenable to an OPEHAT style analysis, since the transfer function varies in a transcen- 
dental (non- algebraic) fashion with the gain value K\ closed form perturbation analysis of the equivalent 
discrete time system is intractable since the matrices A and BC do not commute in general, = 

e eBC = e sc e A < — > A(BC) — (5C)A.) Hence a single “open-loop” OPEHAT transfer function calculated 
with K = 0 is not meaningful to the user. However, a “transcendental root-locus” from u to y may be 
computed by computing the poles and zeros of an equivalent discrete-time plant for various values of k . 

Closed-form analysis of transcendental root-locii is hindered by the property that e (Ac+KB c D c )t cannot 
be computed from e Act and e ( B - c ^ except when A c commutes with B c C c ; this is clearly not the case 
in practice. It may be possible to gain some norm-bounds on the closed-loop eigenvalues of the overall 
discrete-time system in terms of the feedback gain K, but this is an open question. 

4 W-plane analysis 

MARSYAS currently provides z-plane frequency response plots of discrete-time systems. u>-plane analysis 
is often used in practice to allow designers to employ continuous time design techniques to discrete-time 

systems; see, e.g., [ 1 ]. G(z) — ► G(u;) by w = — (^' z ■ s-plane, z-plane, and tu-plane frequency response 

plots are related by z = e* T and w = ^ *.r+} = ^ tanh (^). So as s = jw follows the imaginary axis, 
z = e JwT follows the unit circle (periodic with period 2?r/T), and 

follows the entire imaginary axis each time z = e jwT rotates about the unit circle The inverse bilinear 

transform from w to z is z = ■ — ~ ; hence the Nyquist plots in the z and xv plane will be identical. 

Magnitude and phase information in the 2 and w-planes will be the same except for a “warping” of the 
tu-plane frequency variable v = tan(wT/2). 
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